Parameter identification of a delayed infinite-dimensional heat-exchanger process based on relay feedback and root loci analysis

The focus of this contribution is twofold. The first part aims at the rigorous and complete analysis of pole loci of a simple delayed model, the characteristic function of which is represented by a quasi-polynomial with a non-delay and a delay parameter. The derived spectrum constitutes an infinite set, making it a suitable and simple-enough representative of even high-order process dynamics. The second part intends to apply the simple infinite-dimensional model for relay-based parameter identification of a more complex model of a heating–cooling process with heat exchangers. Processes of this type and construction are widely used in industry. The identification procedure has two substantial steps. The first one adopts the simple model with a low computational effort using the saturated relay that provides a more accurate estimation than the standard on/off test. Then, this result is transformed to the estimation of the initial characteristic equation parameters of the complex infinite-dimensional heat-exchanger model using the exact dominant-pole-loci assignment. The benefit of this technique is that multiple model parameters can be estimated under a single relay test. The second step attempts to estimate the remaining model parameters by various numerical optimization techniques and also to enhance all model parameters via the Autotune Variation Plus relay experiment for comparison. Although the obtained unordinary time and frequency domain responses may yield satisfactory results for control tasks, the identified model parameters may not reflect the actual values of process physical quantities.

2. Saturated relay feedback experiment is performed on the SFOTDM. The detected pole loci are then set to the initial parameters' guess for the HETDM (as a sufficiently accurate mathematical model of the circuit system with heat exchangers), which is enhanced via the LM method, under the single relay-test data. 3. Three scenarios are compared to determine the remaining model parameters and further enhance the already estimated ones via the ATV+ technique and the solution of a nonlinear optimization problem using the LM and NM algorithms. 4. An independent determination of numerator and denominator transfer function coefficients along with the pole loci assignment enables to reduce the number of necessary relay test for some of the scenarios.
The rest of the paper is organized as follows. "Methods and techniques" summarizes theoretical fundamentals of (retarded) TDM spectrum and feedback relay-based experiment, the model parameters' identification using a DF, and the LM and NM methods. "Results" has two fundamental subsections. The first subsection provides the reader with a detailed analysis of the SFOTDM pole loci. The second one presents the HETDM and all steps of determining its parameter values. Namely, the mathematical model of the HETDM is introduced, then the reader is acquainted with the poles assignment, the transfer function numerator estimation using a single relay test, and the complete model parameters estimation via the ATV + technique. In "Discussion", the obtained results are discussed. Finally, "Conclusions" concludes the paper.
The standard notation is used throughout the paper, i.e., C, N, R denote the sets of complex, natural (excluding zero) and real numbers, respectively, R n + expresses the n-dimensional Euclidean space of positive real-valued vectors, Re(s) and Im(s) mean the real and imaginary parts of some s ∈ C , respectively. Superscript T denotes the vector (matrix) transpose.

Methods and techniques
Retarded quasi-polynomial and its spectrum. Let us concisely introduce the Retarded Quasi-Polynomial (RQP) form and its spectrum, i.e., the zero points 9,14,15 . A RQP has the following form where s ∈ C is the Laplace transform variable, q ij ∈ R are non-delay parameters, τ ij ∈ R + with τ i0 = 0 represent delays, and n ∈ N means the RQP order (of derivative).

Definition 1
The RQP spectrum is the set of RQP zeros, i.e.,
2. RQP zeros s k ∈ are isolated and function R n−1 i=0 m i � q ij , τ ij � → s k ∈ C is continuous. 3. For any finite γ ∈ R , the subset � R = {s ∈ � : Res > γ } is finite, while � L = {s ∈ � : Res ≤ γ } is infinite. □ Note that the relation q ij , τ ij → s k is not necessarily smooth; namely, in points where a multiple real root bifurcates into a complex pair.

Definition 2 The RQP spectral abscissa is defined as
Relay-based parameter identification. As introduced above, experimental plant identification using the relay (or another simple nonlinear element) method represents a widely used technique in various engineering and industrial applications. Consider a plant (the model of which is to be identified) under a relay feedback control, as depicted in Fig. 1. In the figure, r(t), e(t), u(t) , and y(t) mean the reference, control error, manipu-  www.nature.com/scientificreports/ lated input and controlled output variables, and G(s) stands for the actual plant (process) dynamics. The choice of r(t) (usually of a constant value) enables to set the operating point.
If the relay parameters are appropriately set, the closed-loop system reaches sustained oscillations of period T osc in a finite time. If the relay element does not cause a phase lag, the corresponding angular frequency ω osc = 2π/T osc is supposed to be close to the so-called ultimate frequency ω u for which ImG jω u = 0 (more precisely ∡G jω u = −π ), j 2 = −1 . However, as a model G m (s) cannot express the true dynamics G(s) exactly, it generally holds that ω osc = ω u . Whenever the relay exposes a phase lag, ω osc < ω u .
By adopting the idea of the DF, one point G jω osc ∈ C can be estimated, i.e., two parameters of G m (s) can be determined. The relay DF, N(·) ∈ C , can be considered as a linear approximation of the relay gain. It is usually derived using a consideration that e(t) is a harmonic signal and u(t) is subject to a truncated Fourier series expansion. Then, for the sustained oscillations, it holds that which enables to estimate parameters of G m (s) . Note that (4) can be graphically interpreted as the intersection of the Nyquist plot of G m (s) with the horizontal line −N −1 (·) . The DF depends on the amplitude A of e(t) oscillations and some other relay setting parameters.
On/off relay test. Let us consider an asymmetrical biased two-level relay. Its static characteristics and the corresponding sustained oscillations (limit cycles) are displayed in Figs. 2 and 3, respectively.
In practice, the setting ε = 0 is suitable when the feedback signal is affected by noise so that the switching relay rate can be reduced. The advantage of the option δ = 0 lies, i.a., in the possibility to estimate the process static gain k = G m (0) = G(0) as (4) N(·)G m jω osc = −1 + 0j ⇔ N(·)G m jω osc = 1 ∡ N(·)G m jω osc = −π  www.nature.com/scientificreports/ for t 0 satisfying that sustained oscillations start for some t < t 0 35 . Purposefully induced asymmetry can also be used to estimate and attenuate the load disturbance 42 . However, it may stop oscillations so that the relay test fails. In addition, model parameter identification with asymmetric relay yields an estimation error of up to 40% in a FOPDT case 50 . Relay with saturation. Estimating the critical point at ω u (or any other ω osc ) does not provide an accurate enough parameter estimation for some processes, e.g., for those with significant time delays. For instance, an error of 23% for FOPDT models was reported 59 .
Model parameters identification can be improved by using saturation relay 35,59 . Its static characteristics and a sketch of the corresponding sustained oscillations (under the assumption of a harmonic output variable) are depicted in Figs. 4 and 5.
The saturation relay does not cause an abrupt step change at e(t) = ±ε , yet it provides a smooth transient around zero. The relay input e(t) is multiplied by k sat resulting in the relay output u(t) up to the limit B = k sat A . The corresponding DF reads Ideally, if the gain k sat is set optimally (i.e., A = A ), input and output signals has the same shape; hence, the DF N sat A, A = k sat is exact. However, in real conditions, u(t) has a shape of the truncated sinusoidal wave with upper and lower limits. Not that the limit case k sat → ∞ yields the ideal relay, i.e. N sat (A, 0) = N(A, 0, 0).
A saturation relay test should follow the standard relay experiment (see the preceding subsection). Once k osc = N(A, ·) is found, then it is set k sat = k min k osc , k min > 1 . Originally, it was suggested to take k min = 1.4 59 . The higher the value is, the closer to the two-level signal u(t) is. Contrariwise, smaller values of k min force u(t) to be closer to sinus-like waves; however, the relay takes a longer time or even can fail to generate sustained oscillation. Figure 3. Sustained oscillations using the asymmetrical biased relay.   38,69,71 introduces an artificial delay τ a > 0 to the serial link between the relay and the process. Every single value of τ a causes the phase lag of ϕ a = ω osc τ a where ω osc means the corresponding angular frequency of sustained oscillations when the delay applies here. Then, the overall phase shift attributed to the process reads Hence, by detecting ω osc and the corresponding amplitude A , another point on the process (model) Nyquist curve can be determined. Obviously, whenever a number of N model parameters are needed to be resolved, then ⌈N/2 − 1⌉ distinct τ a s is required, where ⌈·⌉ means the ceiling function (i.e., the rounding upward to the nearest integer).
The original setting 69 comes from the following idea. The goal is to identify a point located at 45° distance from the negative real axis, i.e., ϕ a = π/4 (under the assumption that ∡N(·) = 0 ). This point is expected to occur at frequency ω osc = 3/5ω u . It eventually yields the following condition and the setting result The disadvantage of this technique is the prolongation of the relay feedback experiment. However, if the initial conditions are sustained oscillations, it lasts a significantly shorter time to restore the oscillations than starting from a constant steady state.
Parameter optimization methods. To solve (2) for given roots s k , (4), and (8), two well-established optimization algorithms are adopted. Their concise description to acquaint the reader follows.   www.nature.com/scientificreports/ J means the Jacobian of f with respect to p , k expresses the iteration step, and > 0 is an adjustable parameter (the so-called damping factor) 93 . Solution (11) and (12) attempts to solve the nonlinear least-squares problem, i.e.,

Levenberg-Marquardt method. Consider a set of nonlinear differentiable functions
The value of (the so-called damping factor) may vary during iterations. One of the framework strategies is to decrease its value as the residual sum on the right-hand side of (13) ( Res p ) decreases, and vice versa. Let us introduce the multiplicative factor κ > 0 as k+1 = k κ . Particular choices of 1 , 1 p , and κ are discussed in "Relay-based parameter identification of heat-exchanger process".
A disadvantage of the LM algorithm is that the solution may converge to a local minimum (as for other Newton-type methods) or it may even diverge (especially, if is set inappropriately).

Nelder-Mead method. Assume an unconstrained optimization problem, first
The idea of the NM method 92 is to iteratively search for the optimal solution by moving a variable-shape simplex in the space of p . The simplex vertices represent the so-called test points. Once the initial simplex 1 S = 1 p 1 , 1 p 2 , . . . , 1 p m+1 is selected, its vertices are re-ordered such that f 1 p i ≤ f 1 p i+1 , i = 1, 2, . . . m (i.e., 1 p 1 represents the best solution estimation) and set k = 1 . Then, the center of the subvector kS = k p 1 , k p 2 , . . . , k p m is computed The worst-valued vertex is reflected through k p c as k p r = k p c + γ r k p c − k p n+1 , γ r > 0 . Then, four scenarios can happen: On condition that f k p e < f k p 1 , set k+1 S = k p 1 , k p 2 , . . . , k p m , k p e , else k+1 S = k p 1 , k p 2 , . . . , k p m , k p r . 3. If f k p m ≤ f k p r < f k p m+1 , the outer contraction is done as k p oc = k p c + γ oc k p r − k p c , 0 < γ oc < 1 .
Whenever f k p oc < f k p r , set k+1 S = k p 1 , k p 2 , . . . , k p m , k p oc , else perform the shrinkage as 4. If f k p r ≥ f k p m+1 , compute the inner contraction k p ic = k p c + γ ic k p m+1 − k p c , 0 < γ ic < 1 . On condition that f k p ic < f k p m+1 , set k+1 S = k p 1 , k p 2 , . . . , k p m , k p ic , else shrink the simplex using (16).
Then k = k + 1 , re-order simplex vertices, and calculate (15), etc. If, however, inequality constraints on a subset p ⊆ p are required, one may use the concept of barrier functions. That is, instead of the objective function f p as in (14), the extended function p is subject to the optimization procedure where β > 0 and f b p > 0 must be sufficiently small as soon as all g j p ≪ 0 ; otherwise, the value of f b p increases considerably until f b p → ∞ as g j p → 0 − .

Results
Root loci analysis of the simple quasi-polynomial. In this subsection, a thorough zero loci analysis of the SFOTDM 12 is provided. The derived results then serve for the pole assignment of the HETDM giving rise to the initial parameters setting of its CE (see "Parameter estimation of the heat-exchanger process model via pole assignment"). The model reads Although pole loci properties of the SFOTDM were studied in the past, according to the authors' best knowledge, a complete image and a thorough exact guide on finding the dominant subset of its spectrum has not been www.nature.com/scientificreports/ provided yet. For instance, Marshal, Gorecki, Walton, and Korytowski 94 studied a generalized characteristic RQP of the SFOTDM with relative parameters ( T = 1, � = ϑ/T, θ = τ/T ), and they determined ranges in which the model is asymptotically stable, aperiodic, and periodic. Moreover, intersections of pole loci trajectories with the imaginary axis for the generalized model were determined. Analogous conditions for which the model is stable, overdamped, critically damped, and underdamped were presented in 95 . Asymptotic behavior of pole loci trajectories in infinity and nearby the imaginary axis were also studied in 96 . Hence, our aim is to analyze the solution (or its rightmost subset) of the CE in C and provide the reader with a simple guide how to compute these pole loci. Result (22) can also be formulate as = e −1 . Let us introduce relative real and imaginary parts of a quasipolynomial root as α = −ϑα , ω = ϑω , respectively. Then the range (23) becomes Lemma 2 means that q SFOTDM (s) has the rightmost double real root at α = 1 for = e −1 . 96 . The double dominant (i.e., rightmost) real root s 1,2 = −1/ϑ (i.e., α = 1 ) becomes a complex conjugate pair for � lim δ→0 + = e −1 + δ . Contrariwise, the double real root becomes a pair of single real roots for � lim δ→0 + = e −1 − δ . □ Theorem 1 q SFOTDM (s) has a real dominant zero in the LHP for = 0, e −1 and a complex conjugate rightmost pair in the LHP for � = e −1 , 0.5π , where the particular root abscissa is within the range (23) (or (24), equivalently). □ Proof From Lemma 1, a negative root abscissa exists only for (23). If ranges from 0 to e −1 , the rightmost real root moves from α = 0 to α = 1 due to Lemma 2 and Lemma 3. From Lemma 2, it is also known that there is the rightmost double real root for = e −1 that bifurcates into a conjugate pair for � = e −1 , 0.5π . Eventually, this pair reaches the imaginary axis again for � = 0.5π as the only (i.e., the rightmost) quasi-polynomial root pair due to Lemma 1. ■

Lemma 3
In the following part of the subsection, dominant (and other) SFOTDM pole loci are investigated.
As an alternative of the proof, one can easily deduce that (28) and (29) are in contradiction. ■ (20) can have only two real solutions (counting multiplicity). □ Proof Lemma 6 implies from Lemma 5 directly due to the root continuity (see Proposition 1). That is, a complex conjugate zeros of q SFOTDM (s) can bifurcate in a pair of distinct real roots only through a multiple pair. Alternatively, distinct real solutions of (20) satisfy (27). Function α → αe −α is unimodal with local and global maximum in α = 1 and the function value 1/e . This point agrees with Lemma 5. Otherwise, the function has two distinct intersections with the constant function ∈ 0, e −1 for α ∈ (0, ∞] . Hence, there is no real solution of (20) for � > e −1 . The situation is illustrated in Fig. 6.  www.nature.com/scientificreports/ (b) If ∈ 0, e −1 , complex conjugate zeros of q SFOTDM (s) are given by (31) and (32), and single real roots are given by the unique solution pair of (c) If = e −1 , complex conjugate zeros of q SFOTDM (s) are given by (31) and (32), and the multiple real root reads s 1,2 = −ϑ −1 (i.e., α = 1). □ Proof Consider item a) first. From Theorem 1 and Lemma 6, there are no real solutions of (20). Complex conjugate ones have to satisfy i.e., both the real and imaginary parts must be equal to zero

Lemma 6 Equation
After some algebraic manipulation, conditions (35) become By expressing �e α from one equation und substituting it into another one, formula (31) is obtained where singularities are to be denied. Further, the latter equation in (36) gives which yields (32) directly. Naturally, only positive right-hand sides of (37) are admissible to get real α values.
We know from Lemmas 5 and 6 and Theorem 1 that the only possible double real root bifurcates into a complex conjugate pair for � > e −1 and there cannot exist another real root of q SFOTDM (s) . Note that only one root from the pair is sufficient to take due to the symmetry.
Assuming item b), a pair of single real roots exists due to Theorem 1. However, it is the only such a pair according to Lemma 6. A single real root must satisfy the first condition in (26), giving rise to (27), the result of which agrees with (33). However, complex conjugate zeros of q SFOTDM (s) must simultaneously exist due to its transcendental manner.
Regarding item c), the existence of the double real root is given by Lemma 2. Besides, there is no other real root of q SFOTDM (s) due to Lemma 5, yet s i ∈ C\R as in (31) and (32) still exist. ■ Theorem 2 does not consider multiple quasi-polynomial roots s i ∈ C\R . The following proposition verifies that such roots can be neglected. (35) and also It is enough to show that a complex conjugate pair of multiplicity 2 does not exist, i.e., n = 1 . We proof a contradiction; hence, let there exists a double root s i ∈ C\R that has to satisfy which gives

Proposition 2 Equation (20) does not admit a multiple pair solution
The latter formula in (39) has two solutions: �e α = 0 or sin ω = 0 . The first one is in the contradiction to the former condition in (39), whereas the second one yields  (39), one gets �e α = 1 , which implies ω = sin (ω) from (36) or (37). That is, the unique solution ω = 0 means that the quasi-polynomial root is real, which gives a contradiction. ■
(Necessity.) Now, we prove by contradiction that ω remains within the limit. Consider that α = (0, 1) . Let exist ω = 0 or ω = π 2 such that the first equation in (41) holds. The limit values are, respectively, which is in contradiction to α = (0, 1). The case ω < 0 can be omitted due to the root loci symmetry in C . Whenever ω > 0.5π , then there exists tan ω < 0 , which implies α < 0 from (31), and we have a contradiction again.
Item c) represents a reformulation of Lemma 2. ■  To conclude this subsection, a graphical procedure to find all the roots of q SFOTDM (s) or the rightmost spectrum of the SFOTDM poles follows. Whenever condition of item b) of Theorem 2 is satisfied, real poles are found as per Fig. 6. Complex conjugate poles are given by (31) and (32) or (41), which can be graphically interpreted as intersections of real-valued functions α 1 (ω) = ω/ tan ω and α 2 (ω) = ln (ω/(� sin ω)) where values α 2 (ω) ∈ C\R determine "forbidden regions" of the function graph, see Relay-based parameter identification of heat-exchanger process. Infinite-dimensional heat-exchanger process model. The HETDM serves as a simulation testbed. The mathematical model arises from heat and mass balance equations that include delays and a thorough analysis of static and dynamic responses of the particular laboratory appliance (see Fig. 8). A concise description of the apparatus follows 84 first. Positions in the figure correspond to the numbers in curly brackets.
The heat fluid circulates in the closed loop flowing through an instantaneous heater {1}, a long insulated coiled pipeline {2}, and a cooler {3}. The power input to the heater (that can be viewed as a solid-liquid flow heat exchanger) is continuously controlled in the pulse-width-modulation sense. Its maximum value is 750 W. The heated fluid temperature on the heater output {4} is only slightly affected when flowing through the 15 m long pipeline; however, the most significant loop delay is caused therein. The outlet temperature of the pipeline is measured by a platinum resistance Pt1000 thermometer {5}. The cooler is constructed as a radiator (i.e., a www.nature.com/scientificreports/ plate-and-fin heat exchanger) that can be considered as an indirect unmixed cross-flow heat exchanger from the process point of view. It is equipped with two cooling fans {6} (one of them is continuously controlled, while another is on/off type). The expansion tank compensates for the impact of the water thermal expansion {7}. The outlet temperature from the cooler is measured by Pt1000 again {8}. Finally, the continuously controllable magnetic drive centrifugal pump {9} serves for fluid circulation. Despite its simplicity, the mathematical formulation of the HETDM and especially its dynamic properties are remarkable due to the model transcendental characteristic equation 84 . As the model is multivariable, the relation between the heater power input u(t) and the cooler outlet heat fluid temperature y(t) is selected as the most interesting input-output pair. Note that both quantities are considered as their deviations from a steady state. The analytically modelled linearized relation reads which is a DDE where b 0 , b 0τ , a 2 , a 1 , a 0 , a 0ϑ ∈ R and τ , τ 0 , ϑ ∈ R + express input/output delays and the state (internal) delay, respectively. The corresponding transfer function is the denominator of which represents the model characteristic RQP (i.e., q HETDM ). In 84 , the following parameter values have been determined by a thorough and complex analysis of static and dynamic data Let us take these data as a benchmark for the significantly more straightforward relay-based experiment. As the values arise from determined physical quantities of the process, they are assumed to be closed to the actual (true) real-life values.

Remark 2
The used Pt1000 thermometers have the guaranteed time constant T 63 of 8 s, i.e., T 90 ≈ 18.4 s. This additional dynamical latency has not been taken into account in analytically-derived model (44), and the true temperature values can be different from the measured ones. However, such negligence does not pose a serious problem with the model. First, plant delays τ , ϑ caused by the thermal fluid transportation have significantly higher values, approx. ≈ 150 s. This means that possible sensor latencies have only a minor effect on the overall dynamics. Second, sensors' latencies do not affect the internal delay of the model itself since they act only in the input/output relation; yet, they are included in the internal delay of the relay-feedback closed loop. If the system is considered linear (indeed, model (44) is a linearized formulation valid in the vicinity of an operating point), a sensor delay can be considered as the additional input/output delay of the model. As input/output delays are not derived analytically but based on measurements, the relay experiment data's evaluation also covers these non-modeled latencies. Therefore, once the model is used for plant control, the plant model and the output signal for the feedback have the same value (in the ideal case).
Simple model parameter estimation using the relay-based experiment. The first step of the identification chain is estimating the SFOTDM parameters, especially those of q SFOTDM (20). It has three substeps. First, the on/off relay with δ > 0 (and ε being sufficiently small), see (5), is used to estimate the static gain k in (19) as per (6). Second, the ideal relay ( δ = 0 ) is applied to get the initial estimation of oscillation data and the input/output delay value. Finally, the saturation relay is used to improve the accuracy of oscillation parameters, which yields the SFOTDM parameters from (4) and (7). All the substeps can be done within a single experiment, saving time since the transition from particular substantial oscillations to others takes less time than setting the oscillations from a constant steady state. Let us use the notation τ → τ s , ϑ → ϑ s for (19) to distinguish the SFOTDM parameters from those of the HETDM (for which no subscript is used). The combination of (4) and (7) can be solved analytically yielding 12 where N · (A, ·) stands for either (5) or (7). However, it is inherently expected that the argument of cos −1 (·) is within the range [−1, 1] . Whenever it does not hold, a numerical solution of the combination of (4) and (7) have to be used instead of (47).
Set B = 100 , δ = 0.05 , and ε = 10 −5 . The relay-test responses are displayed in Fig. 9. The arrows indicate when a particular relay starts to be used.
The eventual data from Fig. 9 are summarized in Table 1. Formula (19) gives k = 3.22 × 10 −2 . The value of τ s can be estimated as the time interval between the switching point of u(t) and the peak time instant of y(t) . Hence, it can be measured that τ s ≈ 136.7 s. Note that k min = 1.4 has been taken for the saturation relay setting, which gives rise to k sat = 185.1 , A = 0.555.
As kN · (A, ·) cos (ω osc τ s ) = −2.72 (for the relay with saturation), (47) cannot be used. Hence, we attempt to apply the NM method to solve the minimization problem www.nature.com/scientificreports/ where G SFOTDM (s) is given (19), N sat A, A and ω osc are taken from Table 1. The barrier function is chosen as f b (T, ϑ s , τ s ) = − x∈{T,ϑ s ,τ s } ln 1 − e −x . The NM control parameters are set to γ r = 1 , γ e = 2 , γ oc = γ ic = γ s = 0.5 . The initial estimation reads (48) [T, ϑ s , τ s ] * = arg min f (T, ϑ s , τ s )  www.nature.com/scientificreports/ The setting of 1 T in (49) arises from the assumption that the inverse of ω osc ≈ ω u is close to the time constant of the delay-free system. The value of 1 ϑ s represents the mid-point of the stability interval (21). Two scenarios for many setting combinations of the initial simplex size and β in (18) are made. First, the fixed value τ s = 136.7 is assumed (i.e., only T, ϑ s are optimized). Second, all three parameters in (48) are to be found.
Among dozens of results (models), six of the most distinguished ones are summarized in Table 2. The results either minimize f (T, ϑ s , τ s ) in the frequency domain or integral absolute error (IAE) and integral time absolute error (ITAE) criteria in the time domain, or represent a trade-off of all the criteria values. Unit step responses (i.e., u = 1 W) are displayed in Fig. 10. The models in Table 2 are also equipped with Roman numerals to label them.
Results IV, V, and VI provide an outstanding cost function value in the frequency domain, yet a worse IAE criterion compared to the remaining results. However, results V and VI give the best ITAE value. Unsatisfactory www.nature.com/scientificreports/ time-domain response of parameters IV can be due to an excellent characterization of the ultimate point G SFOTDM jω u = −N −1 sat A, A but an erroneous estimation of the remaining Nyquist plot. Hence, we do let provide the reader with graphical comparisons of process and model Nyquist plots, see Fig. 11.
Note that the terminal frequency value in Fig. 11 is ω fin = 0.1 rad/s. As can be seen, time and frequency responses for models I, II, and III, and those for models V and VI almost coincide. Figure 11 also indicates that the fixed input/output delay value yields a worse estimation of the ultimate point (i.e., the model crossing point with the negative real axis is quite far from that of the original process characteristics). Moreover, results V and VI seem to give the closest frequency-domain responses to the original Nyquist plot. Hence, Table 3 provides root means squares (RMS) values to measure the error between the process and the models. Two terminal frequency values for RMS computation are chosen, ω fin = 0.1 and ω fin = ω osc = 1.7 × 10 −2 rad/s.  www.nature.com/scientificreports/ Data in Table 3 prove the above-introduced assumption that Models V and VI estimate the process frequency response reasonably on low frequencies (i.e., up to the oscillation frequency); however, they fail for higher ones. It is also worth noting that the estimation of the ultimate frequency is quite accurate since the true one is about ω u = 1.67 × 10 −2 rad/s. Table 2 serve as the initial estimates for HETDM parameters identification. Several scenarios are tested and compared after the HETDM dominant spectrum assignment according to the pole loci of the SFOTDM.

Parameter estimation of the heat-exchanger process model via pole assignment. Now, models (results) in
Initial denominator parameters estimation. The rightmost spectra of SFOTDMs (see Table 2) poles are displayed in Table 4. These loci are computed via the technique presented in "Root loci analysis of the simple quasipolynomial". From Table 4, it can be deduced that the rightmost pole has a decisive impact on the dynamic properties (by comparing the result for SFOTDM I, II, and III), since although the remaining spectra significantly differ, the model time is the same and frequency domain responses are almost identical. Besides, by comparing spectra of models III and VI that are very close to each other, different dynamic properties indicate a high impact of the value of τ s .
The pole assignment problem can be characterized by the set of nonlinear algebraic equations to be solved Since problem (50) includes five unknown parameters to be determined, the unique solution (under the assumption that the Jacobian of the q HETDM has the full rank) requires taking {s 1 , s 2 , s 3 , s 4 , s 5 } . However, it is not always possible. For instance, SFOTDM I has only two significant real poles (the rest are too far in the LHP), and complex conjugate pairs of other spectra cannot be decoupled. Hence, only less number than five of SFOTDM poles are taken in (50). It is also worth noting that whenever s k ∈ C\R , the q HETDM in (50) is split into its real and imaginary parts.
We use the LM method (see "Levenberg-Marquardt method") for solving the pole assignment problem. Let us discuss the ideal initial parameters' selection, 1 p = 1 (a 2 , a 1 , a 0 , a 0ϑ , ϑ) T . Take the CE HETDM and divide both the sides by a 0ϑ Then, it is apparent that the setting yields the CE SFOTDM . However, such a solution is non-feasible and may result in an immature solution of (50). Therefore, several approximate feasible initial settings are eventually chosen.
The multiplicative parameter κ that alternates the damping factor in iteration steps should increase when the residual sum on the right-hand side of (13) increases and vice versa. Some authors suggest setting an asymmetric 97 . Hence, after some numerical experiments, we eventually set κ = 5 when stepping up and κ = 10 when stepping down.
Due to numerous settings of 1 , 1 p , a bunch of possible results is obtained. Selected results (i.e., q HETDM parameter values) and their residual sums Res p opt (see (13)) are provided to the reader in Table 5. Corresponding (50) CE HETDM (s, a 2 , a 1 , a 0 , a 0ϑ , ϑ)| s={s 1 ,s 2 ,...}∈� SFOTDM : q HETDM (s, a 2 , a 1 , a 0 , a 0ϑ , ϑ) = 0 s={s 1 ,s 2 ,...}∈� SFOTDM (51) 1 a 0ϑ s 3 + a 2 a 0ϑ s 2 + a 1 a 0ϑ s + a 0 a 0ϑ + e −ϑs = 0 (52) a 2 = a 0 → 0, a 1 = a 0ϑ T, a 0ϑ → ∞, ϑ = ϑ s www.nature.com/scientificreports/ dominant pole loci of the models are enumerated in Table 6. Note that the QPmR software package 18 is utilized to compute the poles here. It can be deduced that despite diverse results within each of six models (characteristic RQP) families (i.e., I to VI), the obtained dominant spectra are very close to each other. This yields multimodality of the optimization problem. Characteristic RQPs I-a to I-d, II-b, III-a, III-d, III-e, IV-b, V-c, and VI-a to VI-f give almost identical spectra to the original ones (i.e., those of SFOTDMs). These findings correspond to the values of Res p opt in Table 5. It is exciting to observe that even if a subset of two or four poles is prescribed, other one or two complex conjugate pairs coincide with the original spectrum. This feature proves the success of the assignment task and an excellent mapping between q SFOTDM and q HETDM parameters. Now three different scenarios of how to set numerator parameters of the model transfer function (45) or even eventually alter all the transfer function parameters based on relay-experiment data follow.
Numerator parameters estimation using Levenberg-Marquardt method. The first scenario adopts the LM method and utilizes data solely from the single relay test (i.e., without a necessity to perform additional experiments) that should indicate the ultimate data (i.e., the critical point of the Nyquist curve). Once the parameters are determined as in "Initial denominator parameters estimation", they are fixed. Hence, the HETDM transfer function numerator parameters are estimated to comply with conditions (4). The advantage of this scenario lies in fact that all the model parameters are found within a single test.
In more detail, the set of nonlinear algebraic equations to be solved reads where ω osc holds when using the relay with saturation. As only two parameters can be determined by solving (53), two other ones need to be set a priori. Let τ = τ s be fixed as in Table 2. As the static gain k = 3.22 × 10 −2 is known, the following equality is substituted to (53) G o,HETDM (s, ·) = N sat A, A G HETDM (s, ·)  www.nature.com/scientificreports/ Hence, the parameter set in (53) eventually reads p = (b 0τ , τ 0 ) T and b 0 is then calculated using (54). The LM control parameters are set as in "Initial denominator parameters estimation". Table 7 displays the most distinguished results such that two parameters ' sets from each of the six pole spectrum families are selected. The corresponding Res p opt , IAE and ITAE criteria values from unit step responses, and RMS values from Nyquist plots for ω fin = 0.1 and ω fin = 1.7 × 10 −2 are provided in Table 8.
Compared to SFOTDMs (Tables 2, 3), IAE and ITAE criteria values of unit step responses for HETDMs have been enhanced in all six families of models. Contrariwise, frequency-domain error measures have not been improved; the results in Table 8 are very close to those displayed in Table 3.
The unit step responses and Nyquist plots of some selected HETDMs are displayed in Figs. 12 and 13, respectively. The models are selected such they dynamic responses differ significantly. Regarding non-displayed responses, they are very close to some displayed ones. Namely, in the time domain, II-a-1 almost meets III-a-1, yet the latter is faster. II-d-1 and III-c-1 are close to I-a-1. IV-a-1 and IV-b-1 almost coincide, yet IV-b-1 is faster. The same assertion holds for the pair V-c-1 and V-a-1. Finally, VI-e-1 is nearly identical to V-a-1. In the frequency domain, II-a-1 is close to the pair I-a-1 and III-a-1 at the whole frequency range. II-d-1 approaches I-a-1 at low frequencies and I-b-1 at higher ones. III-c-1 is almost identical to III-a-1 for all frequencies. Finally, pairs V-a-1/V-c-1 and VI-c-1/VI-e-1 have frequency responses very close to each other.  Table 6. Dominant pole HETDM spectra.  www.nature.com/scientificreports/ Numerator parameters estimation using Nelder-Mead method from single relay test. Now, let us solve the same task as in the preceding subsection via the NM method. The optimization problem is formulated as follows where the HETDM characteristic RQPs are fixed as in Table 5 and τ = τ s . Again, the value of ω osc is taken from the saturation relay test. Note that the cost function with real and imaginary parts of G o,HETDM (s) is used in (55) rather than that with the amplitude and phase since numerical experiments give better results (in the sense of (4)).
The most outstanding results are provided in Table 9 (two parameters' sets from each of the six pole spectrum families are selected again). Table 10 displays corresponding performance measures and implies that the accuracy of HETDMs in Table 9 is very close (or slightly worse) to models obtained in "Numerator parameters estimation using Levenberg-Marquardt method" (except for models from family I). This means that the models give better accuracy than SFOTDMs in the time domain yet only comparable ones in the frequency domain.
Unit step responses and Nyquist plots of selected models are given in Figs. 14 and 15, respectively. Note again that other characteristics are close to some displayed ones. In the time domain, responses for model families I, II, and III almost coincide when I-a-2 is the fastest and II-c-2 is the slowest. The same assertion holds for families V and VI (VI-b-2 and V-a-2 give the fastest and the slowest response, respectively), while models IV significantly differ from the others. Nyquist plots of III-c-2 and III-f-2 are closest to each other at low frequencies and the (55)   www.nature.com/scientificreports/  Table 9). www.nature.com/scientificreports/ pair II-c-2/III-f-2 at high ones. The pair VI-b-2/VI-d-2 almost coincides at the whole frequency range. Notice that although perfect optimization based on the measured ultimate data (see Table 10) is reached, model IV-a-2 does not provide excellent frequency response. Figure 15 indicates that there must exist a frequency warping. Indeed, the RMS error value of the Nyquist plot for model V-a-2 and ω ∈ 0, 1.7 × 10 −2 is better than that for model VI-b-2. However, the curve for the latter model is closer to the original Nyquist plot than that for the former one. This implies that the model accuracy cannot be judged solely on the shape of characteristics.
Another question is whether optimizing more than three transfer function numerator parameters can improve the HETDMs accuracy.
Numerator/denominator parameters estimation using Nelder-Mead method via Autotune Variation Plus experiment. By substituting (54) into (45), the 8-parameter model is ready to be identified, i.e., q HETDM is not fixed. We do let use the ATV + technique (see "ATV + technique") that dictates the use of three artificial delays yielding the estimation of three additional critical points in the frequency domain. Hence, let τ a,2 = 77.021 s as per (9) and take linear values in the neighborhood of this delay as,τ a,1 = 61.617 , τ a,3 = 92.425 . The common saturated-relayfeedback experiment (see Fig. 16) yields A 2 = 1.460 • C, T osc,2 = 560.03 s , A 1 = 1.328 • C, T osc,1 = 527.23 s , and A 3 = 1.570 • C , T osc,3 = 591.56 s , respectively. Note that k sat = 185.1 , A = 0.555 as in "Simple model parameter estimation using the relay-based experiment".
The NM method is hence used to solve the problem www.nature.com/scientificreports/ where i = 0 is associated with τ a,0 = 0 . Conditional inequalities are incorporated via the barrier function f b (τ 0 , τ , a 2 , a 1 , ϑ) = − x∈{τ 0 ,τ ,a 2 ,a 1 ,ϑ} ln 1 − e −x . If any inequality is broken, the model is not stable or feasible. The remaining parameters, however, can be non-positive. The standard setting γ r = 1 , γ e = 2 , γ oc = γ ic = γ s = 0.5 is adopted while different values of β-see (18)and the initial simplex size are benchmarked. The initial RQP parameter estimates come from Table 5 and let and the value of 1 τ depends on the particular model family, see Table 2, i.e., 1 τ = 136.7, 79.377, 106.869 , or 106.39 s. Results are summarized in Table 11 (again, two eventual HETDMs from each family of models are taken). Corresponding performance measures in the time and frequency domains are provided to the reader in Table 12, and particular dynamic responses for selected models are displayed in Figs. 17 and 18.
Apparently, HETDM parameter identification based on the estimation of four Nyquist plot points results in significantly improved time-domain responses compared to 3-parameter optimizations (see "Numerator parameters estimation using Levenberg-Marquardt method" and "Numerator parameters estimation using Nelder-Mead (56) p * = [b 0 , b 0τ , τ 0 , τ , a 2 , a 1 , a 0 , a 0ϑ , ϑ] * = arg min f (b 0 , b 0τ , τ 0 , τ , a 2 , a 1 , a 0 , a 0ϑ , ϑ) www.nature.com/scientificreports/ method from single relay test"). Regarding the model accuracy in the frequency domain, substantial enhancement is achieved for low frequencies, which, however, does not hold for the whole frequency range. As can be deduced from Tables 8, 10, and 12, the better RMS value for ω fin = 1.7 × 10 −2 is, the worse value for ω fin = 0.1 is and vice versa (except for model family IV). In Fig. 17, nearly nonsmooth step response of HETDM I-a-3 can be observed. An oscillatory response with high-frequency modes of nonnegligible amplitude can be seen for HETDM 4-a-3A.
The remaining (non-displayed) dynamic responses have a similar shape to the displayed ones; however, the plots are not as close as the 3-parameter optimizations. An exception appears for III-f-3, V-c-3, and VI-a-3 where nonminimum-phase-like time responses appear, yet the corresponding Nyquist plots do not prove this feature. The significant step response differences come from diverse input-output delays. Regarding Nyquist Table 11. HETDM transfer function parameters computed via the NM method (artificial delays are used). b0, b 0τ , τ 0 , τ , a 2 , a 1 , a 0 , a 0ϑ   www.nature.com/scientificreports/ plots, the closest curves can be observed for pairs I-a-3/I-b-3 and IV-a-3A/V-b-3 at low frequencies. At higher frequencies, responses differ more significantly, especially those for HETDMs III-c-3, III-f-3, and V-c-3 are far from the remaining ones.
To sum up, most of the models obtained by the solution of task (56) based on the (quadruple) relay-feedback ATV + test gives satisfactory results from the identification point of view.

Discussion
Let us discuss observations made during the entire HETDM identification procedure and also point out some practical issues.
In our experiments, we have supposed that the hysteresis of the on/off relay is negligible, ε ≈ 0 . However, a suitable nonzero value has to be set in practice due to the measurement noise. Such a setting prevents the relay from switching too frequently, which may cause failures. Another important practical issue is the static gain estimation, according to (6). Whenever an asymmetry ( B + � = B − ) is induced, the output of the feedback system also becomes asymmetrical. The problem is that the original output setpoint then shifts, which may cause an erroneous estimation of the static gain. The setpoint shifting can be caused by the feedback nature of a relay experiment, process nonlinearities, and/or disturbances. If someone is unsure about the static gain, the step response test can be made. It is worth noting that disturbances also induce asymmetry of the ideal on/off relay experiment. In such a case, various methods of restoring symmetry can be applied 50 .
As a DF generally represents a linear approximation of a nonlinear element, it is impossible to estimate the critical (or another frequency) point exactly by nature. More precisely, neither the found frequency ω osc nor the loci G m jω osc = −1/N(·) meets the particular values of the actual (measured) process frequency response. This implies that even if the solution of (4) is perfect (see, e.g., HETDM V-a-2 in Table 10), the model does not www.nature.com/scientificreports/ provide sufficient results from the identification point of view. Besides, even if one or more points of the Nyquist curve are estimated well, the remaining course of the plot may vary from the desired loci significantly (see, e.g., HETDM I-a-3 in Fig. 18). As can be seen from Tables 7, 9, and 11, relatively high ratios |b 0 /b 0τ | , |a 0 /a 0ϑ | , |a 0ϑ /a 0 | , (|b 0 | + |b 0τ |)/(b 0 + b 0τ ) often occur. This unpleasant feature yields erroneous steady state computation or numerical instability when simulations (i.e., solving differential equations) due to digital representation of values in computer. Therefore, only some eventual models can be used for control system design and its verification.
Displayed step responses indicate that the initial input/output delay estimation τ = τ s = 136.7 is pretty good. It can be observed that corresponding model families I, II, and III (see Tables 2, 3, 8, and 10) provide slightly better IAE values and significantly lower overall RMS ( ω fin = 0.1 ) compared to models with different values of τ . On the contrary, the better RMS value for low frequencies ( ω fin = 1.7 × 10 −2 ) is, the lower the ITAE value is obtained, which proves the importance of good low-frequency Nyquist plot estimation for the overall timedomain model response.
The eventual SFOTDMs identified using relay feedback experiment were proved to be sufficient for controller design 12 . By matching SFOTDMs dominant pole loci with those of HETDMs and calculating remaining model parameters, performance measures very close to those of SFOTDMs have been obtained in this study. This implies that the eventual HETDMs based on the data from the single relay test (i.e., without artificial delays) can also be used for control tasks. However, the models seem to be insufficient from the identification point of view. Fortunately, the use of the ATV + experiment has brought about much improved models. As many diverse results have been computed with more or less the same cost function values, the identification problem for HETDMS seems to be multimodal task. Therefore, none of the eventual models (see Table 11) approaches the true parameter values (46) given by physical analysis of the thermal process. www.nature.com/scientificreports/ The relay-based experiment can be improved in several ways. For instance, one has to be more careful when setting k sat of the saturated relay (see Fig. 4). The ideal setting should result in sinusoidal-like sustained oscillations of u(t) . However, Figs. 9 and 16 indicate that k sat was set too high. On the other hand, one has to be aware of the necessary existence of sustained oscillations. Another way how to enhance the coefficient value estimation is to capture multiple points of the Nyquist plot, which may yield better matching of process and model curves for a wide frequency range [47][48][49] .

Conclusions
This study should have examined whether it is reasonable to identify parameters of a complex model of a thermal circuit system with internal delays by a parameter identification of a simpler delayed model followed by the models' poles matching. As the identification tool, the standard on/off relay with biased and unbiased feedback test and the relay with saturation have been used. The latter relay should have yielded a more accurate estimation of points on the frequency curve corresponding to sustained oscillations data. Once the simple model is found under a single feedback experiment, its dominant pole loci (of an infinite spectrum) are matched to those of the complex model, giving rise to the characteristic quasi-polynomial coefficients. A simple graphical-based method has been analytically derived to find these loci. The Levenberg-Marquardt method has been applied to solve the pole assignment task. Surprisingly, although only a few poles have been prescribed, some other uncontrolled poles have also been matched. Based on the single-test data, the remaining model parameters have been estimated by the solution of a nonlinear optimization problem (using the Nelder-Mead vs. the Levenberg-Marquardt methods). It has been proved that both the models have had similar time and frequency domain performances. While the eventual models may be sufficient from the control point of view, they fail regarding the accuracy of identified parameters. On the other hand, the proposed procedure enables the estimation of multiple parameters under a single relay test, which is its main benefit. However, we have also performed the ATV + test with artificial delays to get multiple relay-feedback data, which has resulted in much better eventual models.
In the Discussion section, we have touched on some issues that have to be considered in practice and proposed possible further improvements to the proposed concept. Besides, one may apply advanced and more sophisticated optimization approaches to solve the tasks raised in this study. For instance, metaheuristic methods can be benchmarked in the future 98 . In addition, real-life experiments will be made to prove the concept.

Data availability
Data are available by L.P. upon request.